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ABSTRACT 

We describe practical approaches to measuring flexion in observed galaxies. In particular, we look 
at the issues involved in using the Shapelets and HOLICs techniques as means of extracting 2nd order 
lensing information. We also develop an extension of HOLICs to estimate flexion in the presence of 
noise, and with a nearly isotropic PSF. We test both approaches in simple simulated lenses as well as a 
sample of possible background sources from ACS observations of A1689. We find that because noise is 
weighted differently in shapelets and HOLICs approaches, that the correlation between measurements 
of the same object is somewhat diminished, but produce similar scatter due to measurement noise. 

Subject headings: cosmology:observations - galaxies: clusters: general - galaxies: photometry - 
gravitational lensing 



1. INTRODUCTION 
1.1. Motivation 

Flexion has recently been introduced as a means of mea- 
suring small scale variations in weak gravitational lens 
fields (Goldberg & Bacon, 2005; Bacon, Goldberg, Rowe & 
Taylor, 2006, hereafter BGRT). Rather than simply mea- 
suring the ellipticities of arclets, this technique aims to 
measure the "arciness" and "skewness" (collectively re- 
ferred to as "flexion") of a lensed image. Flexion is a 
complementary approach to shear analysis in that it uses 
the odd moments (3rd multipole moments, for example) to 
compute local gradients in a shear field. BGRT have dis- 
cussed how flexion may be used to identify substructure in 
clusters, to normalize the matter power spectrum on sub- 
arcminute scales via "cosmic flexion" (as an analog to cos- 
mic shear), and to estimate the ellipticity of galaxy-galaxy 
lenses. As a practical application, flexion has already been 
used to measure galaxy-galaxy lensing (Goldberg & Bacon, 
2005), and is presently being used in cluster reconstruction 
(Leonard et al., in preparation). 

However, there have been several difficulties in the esti- 
mation of flexion on real objects. First, the flexion inver- 
sion is difficult to describe, contains an enormous number 
of terms, and thus, is rather daunting to code. Secondly, 
there has been little discussion of the explicit effects of 
PSF convolution or deconvolution. Finally, unlike shear, 
there has, until recently, been no simple form to even ap- 
proximate what the "flexion" is. 

The remainder of this paper will thus be a practical 
guide to measuring flexion in real images. We begin, be- 
low, by reminding the reader of the basic terms involved in 
flexion analysis. In § 2, we review shapelet decomposition, 
and discuss some of the issues involved in using shapelets 
to measure flexion. In § 3, we discuss a new, conceptually 
simpler, form of flexion analysis developed by Okura et al. 
(2006), which uses moments, rather than basis functions 
to measure flexion. They call their technique Higher Or- 
der Lensing Image's Characteristics, or HOLICs. We re- 
fine the HOLICs approach somewhat, and develop a KSB 

1 http: / /www. physics. drexel.edu/~goldberg/flexion 



(Kaiser, Squires, & Broadhurst, 1995)-type approach us- 
ing a Gaussian filter to perform an inversion, as well as 
describe a technique for PSF deconvolution. In § 4, we 
discuss comparisons of the two techniques using simulated 
lenses and simulated PSFs. In § 5, we compare shapelets 
and HOLICs inversions on HST images. Finally, in § 6, 
we discuss the implications of this study. 

In Appendix A, we also present the explicit HOLICs in- 
version matrix, so the reader can write his/her own code. 
He/she need not do so, however, as all codes discussed 
herein are available from the flexion webpage. 1 

1.2. Flexion 

What is flexion? Conceptually, flexion represents lo- 
cal variability in the shear field which expresses itself as 
second-order distortions in the coordinate transformation 
between unlensed and lensed images: 



Pi — AijOj + -DijkOjOk, 



(1) 



with 



D ijk = d k Aij , (2) 

where dk is shorthand for d/dxk- Here, A is the normal 
deprojection operator: 



1 - K - 7i -72 
-72 1 - K + 71 



(3) 



and thus, the second term on right in equation (1) repre- 
sents the flexion signal. D may be written as: 



Ai2 



-27i,i — 72,2 —72,1 
—72,1 —72,2 

-72,1 —72,2 
-72,2 271,2 — 72,1 



(4) 



These distortions create asymmetries in a lensed image 
- a skewness and a bending, depending on the values of 
individual coefficients. Irwin & Shmakova (2005;2006) de- 
scribe a similar lensing analysis technique in which the 
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elements of D are referred to as "Catenoids" and "Dis- 
placements." 

BGRT describe an inversion whereby one can estimate 
the individual components, and thus measure two "flex- 
ions" : 

T=d*j (5) 
G = d 1 (6) 
where d is the complex derivative operator: 

d = di + id 2 . (7) 

Figure ?? is reproduced from BGRT and shows the effect 
of a first or second flexion on a circular source. 

An object with first flexion, J 7 , appears skewed, while 
an object with second flexion, Q, appears arced, especially 
if the image has an induced shear as well. The first flexion 
has an m = 1 rotational symmetry, and thus behaves like 
a vector. In particular, it is a direct tracer of the gradient 
of the convergence: 



(8) 



where T\ is the real component and T 2 is the imaginary 
part (as with the second flexion and, as the standard con- 
vention, with shear). 

The second flexion has an m — 3 rotational symmetry, 
though unlike the first flexion, it has no simple physical 
interpretation like that of the first flexion. It is, however, 
roughly proportional to the local derivative of the magni- 
tude of the shear. A more complete discussion of flexion 
formalism can be found in BGRT. 

2. SHAPELETS DECOMPOSITION 

2.1. Review of Shapelets 

Measurement of flexion ultimately requires very accu- 
rate knowledge of the distribution of light in an image. 
The shapelets (Refregier 2003; Refregier & Bacon 2003) 
method of image reconstruction decomposes an image into 
2D Hcrmite polynomial bases: 



(9) 



This technique has a number of very natural advantages. 
In the absence of a PSF, all shapelet coefficients will have 
equal noise. Moreover, the basis set is quite localized (Her- 
mite polynomials have a Gaussian smoothing filter), and 
thus is ideal for modeling galaxies. Furthermore, the gen- 
erating "step-up" and "step-down" operators for the Hcr- 
mite polynomials are simply combinations of the Xi, and 
di operators. 

Refregier (2003) shows that if we decompose a source 
image, /, into shapelet coefficients, the transformation to 
a lensed image may be expressed quite simply as: 



r = {i + nk+ l3 s 3 )j 

where the various lensing operators are: 



(10) 



K = 


1 + I (fit 2 + fit 2 - a\ - fi 2 ) 




Si = 


\ (fil 2 - 4 2 -a\ + fi 2 ) 




s 2 = 


i (a\al - fiia 2 ) , 


(11) 



fit and fi are the normal step-up and step-down operators, 
and the subscript refers to the directional component of 
the coefficient (i.e. 1 for the first or x-component, and 2 for 
the second, or y-component). Note that in the weak field 
limit, these operators indicate that power will be trans- 
ferred between coefficients with indices |An| + |Ato| = 2, 
which preserves symmetry as well as keeping the image 
representation in shapelet space compact. 

In Goldberg & Bacon (2005), similar (albeit more com- 
plicated) transforms were found relating the derivatives of 
shear. We will not reproduce the full second order oper- 
ators here, as they are written in full in the earlier work, 
but we will point out some key features. First, some of the 
elements in the operators have an explicit dependence on 
the (unlcnscd) quadrupolc moments of the light distribu- 
tion. This is due to a relatively subtle effect not present 
in shear analysis. Since the flexion signal is asymmetric, 
the center of brightness in the image plane will no longer 
necessarily correspond to the center of brightness in the 
source plane, and since the shapelet decomposition is per- 
formed around the center of light, we need to correct for 
this. 

Most important, though, is the fact that second order 
lensing terms yield transfer of power between indices with 
|An| + |Am| = 1 or 3. To second order, then, a lensed 
image can be expressed as: 



f' = (l + K K + 1] S,+S^^)f. 



(12) 



Flexion analysis assumes (as does shear analysis) that 
the intrinsic flexion is random, and thus all "odd" (defined 
as n+m) moments are expected to be zero. Thus, from a 
set of shapelet coefficients, a best estimate for the flexion 
signal may be found via x 2 minimization, where: 



x 2 ^ 



Mnimi 



V, 



ll\ 1)1 i //•_> !!!-■ 

(13) 



V nimi n 2 m 2 is the covariance matrix of the shapelet coeffi- 
cients, and [i nm is the "unlensed" estimate of a shapelet 
coefficient. For odd modes, this is zero. For even modes, 
the relative effect of shear is typically much smaller than 
the intrinsic cllipticity of an image, thus it makes sense to 

Set flnm — fnm- 

2.2. Effective Estimation of the Flexion 

Though the form looks quite complicated, conceptually, 
computing the flexion is very straightforward. A simplified 
pipeline may be written as follows: 

1. Generate a catalog of objects and, for each, excise 
an isolated postage stamp. 

2. Compute the shapelet coefficients of the postage 
stamp. 

3. Deconvolve the postage stamp with a known PSF 
kernel. 

4. Compute the transformation matrices associated 
with each of the four flexion operators, solve the % 2 
minimization (equation 13) for jij, and estimate the 
flexion. 
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We discuss each of these steps in turn below. The data 
used for this analysis was taken using HST and the Ad- 
vanced Camera for Surveys, and in the particular context 
of cluster lensing. In this context, the galaxies in which 
we are interested are potentially blended with much larger 
and brighter foreground objects. We discuss the specific 
properties of our data catalog in § 5, but many of the issues 
involved are quite generic. 

2.2.1. Catalog Generation and Postage Stamp Cutout 

The first step in the process, the generation of a cata- 
log and postage stamps seems quite straightforward. For 
some datasets, such as the SDSS (York et al. 2000), the 
data release includes an atlas of pre-cut postage stamps. 
For other applications, such as in relatively shallow galaxy- 
galaxy or cosmic shear /flexion studies, fields will be rela- 
tively uncrowded and thus simple application of widely 
used packages such as SExtractor (Bertin & Arnouts, 
1996) can be used. 

When fields are crowded, however, and contain a wide 
range of brightnesses and sizes, the catalog generation be- 
comes more complicated. It has been noted by Rix et al. 
(2004) that in general, a single set of SExtractor parame- 
ters is insufficient for detection of all the objects of inter- 
est within an image; setting the source detection thresh- 
old too low will result in excessive blending near bright 
objects, whereas a high threshold results in a failure to 
detect fainter sources. Rix et al. describe a two-pass strat- 
egy for object detection and dcblcnding involving an initial 
("cold") pass to identify large, bright objects, followed by 
a lower-threshold ("hot") pass to pick up dimmer objects. 
Their final catalog consists of all the objects detected in 
the cold run, plus any objects detected in the hot run that 
do not lie within the isophotal area of any object detected 
in the first pass. 

This technique works well to prevent spurious deblend- 
ing by SExtractor in images in which there is significant 
substructure. However, when dealing with crowded fields 
(such as clusters of galaxies) the largest problem in cata- 
log generation is excessive blending of sources, particularly 
in the central region. To remedy this, we use a modified 
version of this hot/cold technique. Our method consists 
of a primary SExtractor run to detect only the brightest 
objects. In a lensing field, especially in a lensing cluster, 
these bright objects will tend to be the lenses. Making use 
of the RMS maps generated during this SExtractor run, 
we mask out the bright objects by setting the pixel values 
to background noise, and thus simulate an emptier field. 

We then run SExtractor on the masked image, using 
a much lower detection threshold, to create a catalog 
of background objects. Since shape estimation including 
both flexion and shear have a minimum of 10 degrees of 
freedom, we require at least 10 connected pixels above the 
detection threshold, though in reality, we arc unlikely to 
be able to get a reliable measurement from an image with 
fewer than 15 included pixels. We then discard all objects 
for which reliable shape estimates cannot be found. 

For each remaining object, a postage stamp is gener- 
ated. Ideally, this should identify any neighboring objects 
and mask them out (by setting their pixel values equal to 
background noise). Our postage stamp code also identifies 

2 http:/ /www. astro. caltech.edu/~rjm/shapelets/ 



objects which are blended by using a friends-of-friends al- 
gorithm to find sets of connected pixels that are a certain 
threshold (typically 2-3ct for the stacked images described 
below) above the background. If there is any overlap be- 
tween the object of interest and another object within the 
field of the postage stamp, we consider the source to be 
excessively blended and exclude it from further analysis. 

2.2.2. Shapelet Decomposition 

Shapclcts can be an extremely compact representation 
of an individual image. However, in reality, they are a 
family of basis functions. There is a characteristic scaling 
parameter, /?, which represents the width of the Gaussian 
kernel in the basis function Hermitc polynomials: 



B nm (0) oc exp 



01 +0, 



2/3 2 



(14) 



In principle, while all values of [3 will yield an orthonor- 
mal basis set, some values produce a dramatically faster 
reconstruction in terms of the number of coefficients re- 
quired to reach convergence. Moreover, in reality we don't 
want to reconstruct all details in an image. Structure on 
the individual pixel scale may simply represent noise. 

From a practical perspective, our goal is to optimize se- 
lection of (3, and the maximum coefficient index, n max . 
Rcfregier (2003) suggests the following parameters: 



(3 — \j $ minima v 



n„ 



1 , 



(15) 



where Q m in and 9 max represent the minimum and maxi- 
mum scales of image structure, respectively. 

R. Massey (private communication) has found that 
rather than performing overlap integrals to solve for the 
shapelet coefficients (as was done, for example, in the anal- 
ysis of Goldberg & Bacon 2005), the ideal approach is to 
do a x 2 minimization of the reconstructed image with the 
original postage stamp. This may seem complicated, and 
it is. Fortunately, a shapelets package is available in IDL 
at the shapelets webpage. 2 

For our sample of co-added, background-subtracted, 
HST ACS images of A bcll 1689 , we find that 9 mm = 0.4 
pixels and 9 max = 1.8^/a 2 + b 2 give good shapelet recon- 
structions, where a and b are the semi-major and semi- 
minor axes of the galaxy as measured by SExtractor. How- 
ever, it is important to note that these parameters are 
somewhat dependent on the noise level in the images. 

For a sky-limited sample, we have found that the op- 
timal choice of 9 min scales approximately linearly as the 
ratio of the flux to the RMS sky noise. We have looked 
at this scaling in a sample of galaxies detected with ACS 
(and which we describe in greater detail below), each of 
which was imaged in 4 frames. Prior to stacking of these 
frames, we found that 9 min — 0.75 produced low \ 2 an d 
convergence with small values of n max . After stacking, 
9 min = 0.4 was required. This makes sense, since the nois- 
ier our image, the more prone we might otherwise be to 
fitting complex polynomials to what is, essentially, noise. 

Roughly, the processing time for a decomposition scales 
as Omaxi as ®max determines both the postage stamp size 
and the maximum order of the shapelet decomposition. 
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Due to the high resolution of our images, we encountered 
a number of objects for which n max was so large that the 
decomposition time became prohibitive. We opted to re- 
grid images with n max > 50 into larger pixels by taking 
the mean of the pixel values in square bins, the size of 
which is determined by binsize = n max /50. This number 
was rounded up for objects with 50 < n max < 75. 

A flexion measurement is then carried out using the \ 2 
minimization technique described previously. However, we 
have found that truncating the shapelet series prior to the 
flexion measurement yields a more accurate and robust 
measure of the flexion than using the full series. Excluding 
the higher order shapelet modes avoids contamination of 
the flexion signal by small scale substructure and by noise 
(particularly in dimmer objects). We exclude all shapelet 
modes with n > n max /5 in our flexion measurement. This 
effectively increases 9 m i n to 2 pixels, without affecting the 
accuracy of the reconstruction. 

2.2.3. PSF Deconvolution 

One of the complications in measuring properties of 
lensed images is that, in practice, they are convolved with 
a PSF: 

r d 2 6'P(6-e')f {0 \6) . 



f(0) 



(16) 



In principle, the PSF can be estimated through mea- 
surement of stars, but in deep, small-field, high galactic- 
latitude observations, stars may be scarce, and thus PSF 
estimation may rely partly on numerical analysis of the 
instrument (e.g. the Tiny Tim algorithm, Krist, J. 1993). 

In reality, though, this should rarely be an issue. Es- 
timations of the PSF flexion from Tiny Tim yield values 
of (JaF.psf — 7 x 10 -4 . This represents the maximum in- 
duced flexion which can arise from convolution with the 
PSF, and is still several orders of magnitude lower than 
the scatter in intrinsic flexion of galaxies. 

We are not surprised by this since, for example, PSF 
distortions arising from variable sizes in chips is likely to 
scale as the variation in PSF ellipticity. In ACS, chip 
distortions produce ellipticities of order 1%, and vary on 
scales of 100's of pixels, producing an induced flexion of 
~ 10~ 4 pixr 1 . From the ground, the atmospheric distor- 
tions are expected, on average, to be even more isotropic. 

There is another reason to suppose that PSF flexion con- 
tributions will be unimportant. In shear measurements, 
the PSF ellipticity typically varies smoothly and somewhat 
symmetrically around the center of a field, mimicking (or 
partially reversing) the overall behavior of the expected 
shear field. Since flexion probes smaller scale effects, the 
induced flexion by the PSF will, on average, cancel out. 

This is not to say that we cannot deal with PSF flexion 
inversion. Refregier (2003) describes an explicit deconvo- 
lution algorithm (see also Refregier and Bacon 2003, and 
references therein). In shapelet space, equation (16) can 
be re-written as: 



In 



E °> 



nmn'm / n // m" 



P AO) 

r n'm' J n "rr 



(17) 



n' m'n"m" 



Where C nmn > 

m'n"m" is the 2-dimcnsional convolution 

tensor: 



and a, (3 and 7 are the characteristic scales of f^°\ P and 
/, respectively. We may then define a PSF convolution 
matrix as: 



I'm' — ^ ^ 



; n in' m' n" m" Fn" m" 



(19) 



If only low order terms in the convolution matrix are in- 
cluded, it may be inverted to perform a deconvolution via: 



nmn'm' fn'', 



(20) 



This provides a good estimate of the low order coefficients, 
but high order information is lost. An alternative inversion 
scheme involves fitting the observed galaxy coefficients us- 
ing a x 2 minimization scheme. Refregier and Bacon (2003) 
note that the \ 2 scheme may be more robust numerically, 
and can take full account of variations in the noise charac- 
teristics across an image (although it is strictly only valid 
in the case of Gaussian noise). It is this scheme that is 
implemented in the shapelcts IDL software. 

2.2.4. Flexion Inversion 

If the shapelet coefficients are statistically independent 
(as they will be in the absence of an explicit PSF decon- 
volution), formal inversion of the flexion operator is quite 
straightforward. Under these circumstances, we also have 
the benefit that the measurement error for each moment 
is identical (see Refregier 2003 for discussion). 

Noting that, in most galaxies, the coefficients corre- 
sponding to the n + m =even moments will be much larger 
than the odd moments (and, indeed, upon random rota- 
tions, the latter will necessarily average to zero) we can 
dramatically simplify equation (13). First, we define the 
susceptibility of each odd moment as: 



0(2) f 

J ij Jnm 



where f nm represents all of the "even" coefficients, and 
n'm' represents all of the odd coefficients. Thus, we wish 
to solve for the relation: 



E 



(fn' 



1i,j&f n < m ',ij) 2 = min. 



(21) 



where the first term is taken directly from measurement. 
Taking the derivatives and rearranging, we find: 



m' ,ij A/, 



n'm' ,i'j' j 



(22) 



which can readily be inverted to solve for ~fi t j. 

In practice, however, there are a number of issues which 
must be considered. First, if the PSF or pixel scale are rel- 
atively large compared to the minimum resolution scale of 
an image then many of the high-order moments returned 
by shapelets decomposition will, in fact, not have any in- 
formation. Thus, the above inversion will yield a system- 
atic underestimate of the true image flexion. Above, we 
describe a truncation which minimizes this effect. 



C n mn'm'n"m' 



While the flexion inversion is, at its core, linear alge- 
(18) flexion along with examples on the flexion webpage. 
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3. HOLICS ANALYSIS 

3.1. Higher Order Moments 

Okura et al. (2006) recently related flexion directly to 
the 3rd moments of observed images. This is a signifi- 
cant extension of flexion, and very much along the lines of 
Goldberg & Natarajan's (2002) original work which talked 
about "arciness" in terms of the measured octopole mo- 
ments. Throughout our discussion, we will use the nota- 
tion: 

r {0 i -el)(e j -e])f(e)d 2 6 (23) 



Qij — 



4/< 



to refer, in this case, to the unweighted quadrupole mo- 
ments, with all higher moments being defined by exact 
analogy. In this context, F refers to the unweighted inte- 
grated flux. 

They define the complex terms: 

c ^ (Qui + Q122) + t(gii2 + Q222) (24) 



and 



s = (Qui - 3Q122) + i(3Qn2 - Q222) ^ 



where 

£ = Qnii + 2Q 1122 + Q2222 • (26) 
These terms are collectively referred to as HOLICs. 

If a galaxy is otherwise perfectly circular (i.e. no ellip- 
ticity), and in the absence of noise, then the HOLICs may 
be directly related to estimators of the flexion (subject to 
an unknown bias of 1 — k). Namely: 



T ■ 



9£-6(Q? 1 + Q2 2) 



0^ 
3 



(27) 



(28) 



where the latter term in the denominator of T docs not 
appear in the Okura et al. analysis. Bacon and Goldberg 
(2005) show that a flexion induces a shift in the centroid 
proportional to the quadrupole moments. In order to cor- 
rectly invert the HOLICs, this term needs to be incorpo- 
rated explicitly. The simplicity of the extra term results 
from an approximation of near circularity. 

The beauty of this approach is that it gives us a very 
intuitive feel for what flexion means in an observational 
way. We thus introduce the term "skewness" to the in- 
trinsic properties of a galaxy as measured from equa- 
tion (27) whether or not the galaxy is otherwise circular, 
and whether or not it is lensed. The skewness may be 
thought of as the intrinsic property, much as the "elliptic- 
ity" is the intrinsic property related to the "shear." Like- 
wise, the intrinsic property associated with equation (28) 
will be referred to as the "arciness." 

In reality, however, equations (27) and (28) are not suf- 
ficient to perform a flexion estimate even if a galaxy has 
an ellipticity of only a few percent. Okura et al provide 
a general relationship between estimators for flexion and 
HOLICs, though the relation is best expressed in matrix 
form: 



M 



( k \ 








Qi 




V Q2 J 





(29) 



where M. is a 4 x 4 matrix consisting of elements propor- 
tional to sums of Qijki and QijQki, the former of which 
can be found by explicitly expanding the expressions in 
Okura et al., and the latter of which is again derived from 
the shift in the centroid. For the convenience of the reader, 
we write out the explicit form of M. in Appendix A. 

It may be seen by examining the elements of M. why 
this inversion must be done explicitly for even mildly el- 
liptical sources. For fully circular sources, it may be seen 
by inspection that M. is diagonal. However, when a source 
has an ellipticity even as small as 10%, it can be shown 
that |Mn| ~ |Mi 2 |, and thus equations (27) and (28) are 
no longer even approximately correct. 

3.2. Gaussian Weighting with HOLICs 

The application of the HOLICs technique would be triv- 
ial if there were no measurement noise. In the presence of 
noise, and especially, when the sky dominates, measure- 
ment of unweighted moments is inherently quite noisy. In 
a case where we are measuring the 3rd and 4th moments, 
it is even more so. 

Kaiser, Squires & Broadhurst (1995; see also a nice re- 
view by Bartelmann & Schneider, 2001) developed perhaps 
the most comprehensive approach to dealing with the sec- 
ond moments (the ellipticity) with noisy observing, and 
with a (potentially anisotropic) PSF. 

Our approach is similar. We have only worked with a 
Gaussian window thus far, but the approach is general- 
izable for any circularly symmetric weighting. We thus 
define a window function: 



W{9) 



1 



exp 



(30) 



where the origin is taken to be the center of light, and 
the integral is normalized to unity. Further, we define the 
weighted moments as, for example: 

On = j J (0i - T 1 ) 2 .f(e)w(e)d 2 e . (3i) 

We can thus redefine all HOLICs and moments similarly. 
We have found through experimentation (see below) that 
for a sky noise limited source, a reasonable value of aw is 
1.5 times the half-light radius. 

If we were to simply replace all elements in M, (, and 
S from equation (29) with their weighted counterparts we 
would not get an unbiased estimate of the flexion. There 
are two corrections. One has to do with the fact that 
centroid shift will differ from the unweighted case to the 
weighted case. Consider an extreme scenario in which the 
window width is arbitrarily small and in which the un- 
lensed image was circularly symmetric with a peak at the 
center. In that case, the centroid will essentially remain 
at the center (peak brightness) even if the unweighted mo- 
ments shift. 

Thus, compared to the unweighted moments, the cen- 
troid will shift: 



Mi = 



a 



(32) 
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where we have used the explicit fact that for a Gaussian: 



Making the substitution, 



dW{9) 



■-S-W(0) , 



(33) 



' w 



The other correction has to do with the fact that though 
lensing preserves surface brightness, it does not preserve 
total flux. This is normally related by the Jacobian of the 
coordinate transformation. However, when considering a 
window function, we need consider that transformation 
explicitly: 

dp 



W{/3)d 2 /3 



de 



W(/3)d 2 



(34) 



as used by Okura et al. (2006), and where we have simply 
multiplied both sides by the factor W((3). In this context, 
f3 refers to the image coordinate in the source plane. Ig- 
noring the terms proportional to shear (which cannot be 
directly addressed by this method at any rate), we have 
the approximate relation: 



1 



dW 



W(f3)~W(6) + -D ijk e i 6 j — 



(35) 



or, as we have already asserted: 



W{(5) ~ W{6) - l D ijk e -^W{6) . (36) 
1 a w 

Note that the latter term contains an odd number of posi- 
tion elements, and thus, coupling to the generating equa- 
tions for ( and 8 produces contributions of 6 th moments 
in M: 



AH ——In Qijklmn 



'w 



which, in turn, must be corrected for. 
We may thus say that: 

M = M{Q ij ,...) + AM , 



(37) 



(38) 



where the latter expression can also be found in Appendix 
A. 

3.3. PSF Correction in HOLICs 

As with our discussion of shapelets, above, we must 
also consider PSF deconvolution in our HOLICs pipeline. 
We define the PSF function in equation (16), and all un- 
weighted moments of the PSF are denoted by Py, etc. In 
principle, because of the higher signal-noise of the PSF, 
the unweighted moments are easier to estimate than the 
moments of the detected image. While we argued, above, 
that the induced flexion from a PSF is likely to be small, 
it is still the with shear, that the PSF will reduce 

the measured flexion. Let's first consider the case in which 
we were able to measure the unweighted moments of both 
the PSF and the observed image. It is straightforward to 
show that: 



Qa = Q 



(0) 
ij 



Pi; 



Thus may be computed via the relation: 



Q 



9 i 9 J f m (9')P(9-9')d 2 9d 2 9 / 



(39) 



(40) 



(41) 



yields 



(42) 

It is straightforward to show that this yields equation (39). 
Similarly, it may be shown that: 



However, 



Qak = Q<gL + Piik ■ 



Quu = Qifii + Pun + 6Qi°Mi 



(43) 



(44) 



with a similar expression for Q2222, and 



Q1122 = O1122 + + Qii ) P 2 2 + <^°Mi (45) 

provided we assume the PSF is nearly circular. 

If we further look only at nearly circular sources, then we 
may estimate the flexion using the forms in equations (27) 
and (28). Again, assuming unweighted moments, and zero 
PSF and intrinsic flexion we find: 



T = T 



9gW-6(Q<° )2 + Q^r) 
^ - 6(Q 2 n + Q 2 22 ) 



,(0)2,, 



(46) 



Where T% is an unbiased estimate of the flexion, and Ti is 
the estimated flexion if one does not include the correction 
for the PSF. However, the normalization constant may be 
estimated directly from combinations of the PSF 2nd and 
4th moments, and the unweighted moments of the image. 
Since this term represents something like the overall ra- 
dial profile of the source, the unweighted moments can be 
estimated even under noisy conditions. 

Similarly, the second flexion may be estimated as: 



Gi = Q; 



£(0) 

T 



(47) 



Though we have derived these relations for a nearly cir- 
cular source, we have found they provide a good correction 
even when the PSF and intrinsic image size are compara- 
ble, and when ellipticities for the source image are e ~ 0.2. 

4. SIMULATED LENSING 

Which approach is better, shapelets or HOLICs? From 
a signal perspective, the shapelets technique is better. It 
is designed to provide optimal weighting and return opti- 
mal signal-noise. Moreover, as described above, inversion 
of the PSF is a straightforward and well-designed process. 
In the absence of noise, the two techniques produce very 
similar results. 

On the other hand, the HOLICs technique has several 
practical advantages, especially for large surveys. For one, 
the HOLICs code is typically much faster than shapelets. 
For an N pixel image, the HOLICs technique is an O(N) 
calculation, whereas the shapelets is 0(N 2 ). Addition- 
ally, some values of (3 produce very bad reconstructions, 
and hence, minimization of x 2 can be time-consuming and 
may not converge to a minimum. 
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As a simple test, we created images with brightness pro- 
files of: 

J(r) oc exp[-(r/r )™] (48) 

and though we found similar results for a reasonable range 
of exponents, the results presented below are for n = 1.5. 
We have used a constant source ellipticity, typical of those 
observed in the field, e = 0.2, and had measurement errors 
which were dominated by sky brightness. In each case, we 
had no intrinsic arciness or skewness (that is, the flexion 
of the unlcnsed objects were zero), since our aim was to 
measure the response of each of the estimators to lensing. 

We then artificially lensed each of our simulated im- 
ages, added sky noise, and measured the flexion using both 
the HOLICs and shapelets techniques. The noise is fixed 
throughout this discussion, as is the strength of the flexion 
signal. It is clear, however, that all relevant signal-noise 
values will scale linearly with the strength of the lensing 
signal and inversely with sky noise. 

4.1. Optimzing the HOLICs Scale Factor 

Our first questions is, what is the optimal value of Cw, 
such that: 

<?W = C W X Thai f -light ? (49) 

Ideally, we would like an unbiased estimator of the flexion 
which also has very little scatter. It is clear that the larger 
the value of Cw, the larger the scatter will be (in general), 
since we will be measuring more and more of the noisy sky. 
However, the smaller the Cw, the less accurate will be our 
measure of the real shape of the galaxy. Figure ?? bears 
this out. There is an optimal value of Cw around 1.5, 
which reflects a balance between minimizing measurement 
errors as well as any measurement bias inherent in the 
technique. 

With shapelets, we find a systematic underestimate of 
11% in the first flexion, and an overestimate of 12% in the 
second flexion. We find a scatter of about 12% in both. 
This is very similar in magnitude to the results found by 
an "optimal" HOLICs analysis. 

4.2. Correlation of HOLICs and Shapelets Measurement 

Error 

Since both HOLICs and shapelets give similar measure- 
ment errors at fixed sky noise, it is worth considering 
whether we expect measurement errors between the two 
techniques to be correlated. Even in these idealized cir- 
cumstances, uncorrelated errors would mean that there 
is significant information in the images which is not being 
used. In Fig. ??, we show the correlation in uncertainty be- 
tween our Cw = 1-5 HOLICs estimates, and our shapelets 
estimates. 

For the first flexion, in particular, the correlation is quite 
high, with a Pearson correlation coefficient of 0.86. The 
correlation in measurements of the second flexion is much 
lower, with pq — 0.23. Why don't they have perfect corre- 
lation? The two techniques weight various components of 
the signal (and thus, the noise) differently, and therefore 
have a slightly different response to the noise. 

This general trend is borne out with observed objects as 
well, in which we will see much higher correlation between 

3 http://terapix.iap.fr 



measurements of the first flexion than the second flexion 
between the two techniques. 

4.3. PSF Deconvolution 

Finally, we can simulate PSF deconvolution. Using a 
Gaussian PSF with a characteristic size somewhat larger 
than the intrinsic image (the correction factor described 
in equation 46 is 2.7), we distorted and then recovered the 
flexion estimates from images of increasing intrinsic ellip- 
ticity. This analysis is done in the absense of sky noise, 
and thus any errors in shape recovery represent a system- 
atic effect. We show the fractional errors in measurement 
of the first and second flexion in Figure ??. Since it is pos- 
sible to estimate the systematic error for a combination of 
measured shear and PSF shape, it is advisable to those 
wishing to make high-precision flexion measurements to 
take this empirical correction into account. 

We find that, despite the fact that the PSF correction 
is based on an assumption of circularity, it continues to 
produce a good result even if the image has an intrinsic 
ellipticity as high as 0.3. 

5. MEASUREMENT OF FLEXION ON HST IMAGES 
5.1. Sample Selection and Pipeline 

We also compare the two approaches to flexion inver- 
sion on real objects. Our data consists of 4 HST ACS 
cosmic-ray rejected (CRJ) images of Abell 1689 using the 
F625W WFC filter (hereafter "R-band"). Each image was 
taken by H. Ford during HST Cycle 11, and has an ex- 
posure time between 2300-2400 seconds. The observations 
are described in detail in Broadhurst ct al. (2005). 

Using the SWarp software package 3 , these four images 
were co-added to create a single "full" R-band image. We 
also generated 2 independent "split" images for compar- 
ison purposes by combining only two of the original im- 
ages. The images are background-subtracted, aligned and 
re-sampled, then projected into subsections of the output 
frame using a gnomonic (or tangential) projection, and 
combined using median pixel values. 

Each image undergoes a primary SExtractor run de- 
signed to detect only the foreground objects (cluster mem- 
bers and known stars). This detection is carried out us- 
ing the cross-correlation utility in SExtractor, which al- 
lows us to specify the locations of the foreground objects. 
Our foreground object catalog was generated using a com- 
bination of spectroscopically confirmed cluster members 
(Due et al. 2002), and identification by eye of foreground 
objects that were later confirmed as such by use of the 
NASA/IPAC Extragalactic Database (NED), as well as 
clearly identifiable stars in the field. These objects are 
then masked out as described previously, and a second 
SExtractor run carried out. 

A catalog of objects is then generated, using only those 
objects that were detected in both of the split R-band im- 
ages. We measure the flexion in our catalog objects using 
both the truncated shapelets method (described above) 
and the HOLICs approach, and then compare the mea- 
surements by computing Pearson correlation coefficients 
between the different estimates in the full image. We also 
compute correlation coefficients between measurements 
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taken using the same technique in the two split images. 
This gives us an estimate of the robustness of the mea- 
surement technique. 

When computing the correlation coefficients, we include 
only objects with a > 3 pixels, and consider only the 
brightest half of our catalog objects. In order to exclude 
extreme or erroneous measurements, we require (a|jF|) < 
0.2 and (a\G\) < 0.5. 

5.2. Results 

Figure ?? shows a comparison of the HOLICs and 
shapelets estimates of flexion in the full image. Both T 
and Q have a positive correlation, with a Pearson correla- 
tion coefficient of 0.17 for F and 0.12 for Q. Additionally, 
both methods yield similar standard deviations for both 
first and second flexion. 

This is what we expected from our simulated results 
above. Clearly, if flexion represents any real signal, the two 
techniques should be correlated, and, as we showed in our 
simulated results, the correlation in first flexion is higher 
than in second flexion. But the correlation in our mea- 
sured results is lower than in the simulated ones. Why? 
In part, this is due to a relatively noisy field. We've found 
that selections on brighter magnitudes and larger objects 
improves the correlation somewhat. In part, however, this 
is due to what we mean by "flexion." Recall that the 
shapelets and HOLICs analysis of flexion involve weighting 
different modes in different ways. Real, unlensed, galaxies 
will have odd modes which are not necessarily correlated 
in a simple or obvious way. Lensing, of course, produces 
a significant correlation, and thus, a population of signifi- 
cantly lensed objects (for which the majority of the flexion 
is due to lensing) would be expected to have a more corre- 
lated flexion. This is similar to the case with weak shear 
analysis in that the S/N from a typical object is usually 
less than 1. 

We can test this hypothesis directly by comparing the 
measurements in the split images and estimating the flex- 
ion in both using the same technique. Any discrepancies 
between the two ought to be the result of photon noise 
rather than intrinsic complexity in the structure of the 
3rd moments. 

Figure ?? shows a comparison of the HOLICs measure- 
ments made on each of the split images. These measure- 
ments are well correlated: the Pearson correlation coeffi- 
cient here is 0.37 for T and 0.23 for Q. In Figure ??, we 
see a comparison of the shapelets measurements in these 
images, which appear to be more strongly correlated (par- 
ticularly for T). The Pearson correlation coefficients here 
are 0.58 for T and 0.18 for Q. 



As motivated above, most of the "noise" in our mea- 
surements comes from the intrinsic distribution of flexion 
within our sample. Indeed, using the HOLICs approach, 
we find: 

V a \F\ + Noise = 0-05 (50) 
&a\G\+Noise = 0-08 (51) 

The distribution function may be seen in Fig. ??. Note 
that this result includes noise. However, we may estimate 



the relative effect of photon noise on this scatter by using 
correlation between frames. That is: 

a a\F\ = y/P&a\F\+Noise (52) 

And thus, we find that our best estimate of the intrinsic 
scatter in first flexion is: 

<T a \F\ = 0.03 (53) 

(as found in Goldberg & Bacon 2005), and 

<Ja\G\ = 0-04 . (54) 

The combination a\T\ represents a dimensionless term, 
and thus is independent of distance. 

It should also be noted that since these measurements 
are taken within a cluster, the signal is included as well, 
and one might question whether it is reasonable to esti- 
mate the intrinsic variability of flexion from lensed images. 
The intrinsic scatter in flexion was originally measured in 
Goldberg & Bacon (2005), and we merely confirm the re- 
sult here. However, this is a reasonable thing to do, as 
flexion drops off much more quickly than shear, and thus, 
even within a rich cluster, the flexion signal is dominated 
by individual galaxies. Even at a separation of 1", the 
flexion from even a very massive 300fcm/s galaxy on an 
a = 0.4" source is about 0.05, approximately the level of 
the intrinsic flexion. Such separations are relatively rare, 
however. 

6. DISCUSSION 

We have endeavored to present a detailed guide to mea- 
suring flexion in real observations, with a focus on space- 
based imaging. In the process, we have taken a look at 
two different approaches to measuring flexion: shapelets 
and HOLICs, with an eye toward which approach is "bet- 
ter." From an idealized perspective of maximal signal- 
noise, the answer is simple. Shapelets produces a mode- 
by-modc comparison which optimally averages to produce 
a unique estimate of flexion. However, this result is com- 
plicated somewhat in two limits: blending, which affects 
larger objects, and PSF convolution, which affects smaller 
ones. 

When images are blended, it is clear that we benefit by 
giving extra weighting to those pixels near the center of the 
the object. In that sense, HOLICs can be said to produce 
more robust results. Likewise, despite an explicit PSF de- 
convolution algorithm, applying the flexion inversion using 
shapelets results in inclusion of small scale power which 
has been blended away through the atmosphere or instru- 
ment. We have discussed, above, how this might be alle- 
viated by only using relatively low order modes from the 
reconstruction in the estimate of flexion. However, doing 
so comes at the expense of some (but by no means all), 
of the signal- noise advantage from shapelets. Indeed, even 
using a relatively truncated form of the shapelets analy- 
sis still produced greater correlation between independent 
images of the same objects, and thus, cleaner estimates of 
the flexion. 

However, one complication in the shapelets analysis is 
producing a good shapelets decomposition in the first 
place. While R. Massey's shapelet code comes with an 
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optimization routine to find the best fit scaling parame- 
ter, (3, the shapelet decomposition runs several orders of 
magnitude slower than HOLICs. For very large lensing 
fields, this may prove a significant limitation, and thus, 
HOLICs provides a fast, physically motivated, reasonably 
reliable alternative. 
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APPENDIX 

EXPANDED COEFFICIENTS FOR HOLICS ANALYSIS OF FLEXION 
In equation (29), we state that the flexion may be solved via inversion of the relation: 

y = Mx 



(Al) 



where a; is a vector of the desired flexion estimators, and y is the measure of the 3rd order HOLICs. Here, we show the 
explicit form of M. 



Mn 

Mi 2 
M13 
M14 

M 21 
M 22 
M 23 
M 2 a 
M 3 i 
M 32 

M 33 
M 34 
M4.1 



M, 



(9 + 8j?i 



33Qfi + 14QnQ 22 + Q 2 22 + 20Q 2 



1 

32Qi 2 Q 22 + 32(5iiQi 2 



12 



2r? : 



4£ 



l/ n , \ \ SQii ~ 2Q11Q22 — Q22 ~ 4Qi 2 
4(2* + Ai) ^ 

1/r, \ \ 2QllQl2 

-(2t? 2 + A 2 ) 



= 2772- 



32Qi 2 Q 2 2 + 32QiiQi2 
4£ 



42 



1, s - q\ On + 14Qiig 22 + 20Q 2 2 + 33Qj 2 
-(-8r?i +9) 

l/o \ \ -2Qi 2 Q 22 
-(-2r? 2 + A 2 ) — 

1, . x (Qfi+4Qf 2 + QiiQ 22 -3Qj 2 ) 
-(2^! - Ai) 

l,, n ^ 7i v 3(UQfi - 10QnQ 22 - QI2 - 20Q 2 2 ) 
-(11% + 7Ai) — 

1, in ^ 7 ,^ 3(8QnQi2-32Q 12 Q 22 ) 
^(-lOifc + 7A 2 ) — 

3 3(~2QnQ 22 + Q 2 U + Q 2 22 + 4Q 2 12 ) 

4 4£ 


l,, n ^ 7X v 3(32QnQ 12 - 8Qi 2 Q 22 ) 
-(10r? 2 + 7A 2 ) 

l,, n 7i v HQii + 20Q 2 2 + IOQ11Q22 - llgia) 
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M 43 = 



M 44 = (A2) 

where, as defined in Okura et al. (2006), we use: 

(Q1111 - Q2222) + 2z(Qni2 + Q1222) , a n 

V = £ (A3) 

^ _ {Qnn - 6Q1122 + Q2222) + 4^(Qni2 - Q1222) 

Note that i] = and A = for all circularly symmetric distributions and even those with no cllipticity but with flexion. 

If we apply a Gaussian weighting with width, aw, to our moment measurements, then M. should be computed using 
the weighted moments. In addition, the following terms must be added: 

A „, — 3Q1IIIH — 6QllH22 — 3Qll2222 + 3Q22Qll22 + 9(5llQllll + 6Ql2QlH2 + 9QllQll22 + 6Ql2Ql222 — 3Q22QHH 

AMn — — -s 

— 3Qlllll2 — 6Q111222 — 3Qi22222 + 3Q22Qlll2 + 9<5nQl222 + 3Q22Ql222 + 6Q12Q1122 + 9QllQlll2 + 6Q12Q2222 



AM 21 
AM 22 
AM 23 
AM 24 
AM 31 
AM 32 
AAf 33 
AM34 
AM 41 
AM 42 
AM43 
AM44 



— Qllllll + 2Qnii22 + 3Qn2222 — 3Q22Qll22 + 3(5iiQnii + 2Qi2<5lll2 — 9QnQii22 — 6Q12Q1222 + Q22Q1IH 

— 3Qlllll2 — 2Qiii 2 22 + Ql22222 + 6Q12Q1122 + ^QllQllVl — 3QllQl222 + 3Q22Q1II2 _ Q22Q1222 _ 2Qi2Q 2 222 
~3Qimi2 — 6Q111222 — 3Qi22222 + 6Q12Q1122 + 3QiiQni2 + 9Q 2 2Qlll2 + 3QnQi222 + 9Q 2 2Ql222 + 6Q12Q1IH 

— 3<5nii22 — 6Qn2222 — 3Q222222 + 6Q12Q1112 + 3QiiQ2222 + 6Q12Q1222 + 9Q22Qll22 + 3Qii<5ii22 + 9Q22Q 2 222 

-Q111112 + 2Q111222 + 6Q122222 - 6Q12Q1122 + Q11Q1112 + 3Q22Q1112 - 3Q11Q1222 - 9Q22Q1222 + 2Q12Q1111 

— 3Qllll22 — 2Qii2222 + Q 2 22222 + 9Q22<9ll22 + 3(5iiQn22 — QllQ2222 + 6Q12Q1112 — 2Qi 2 (5l222 — 3Q22<92222 

jQmm + 6Q111122 + 9Q112222 - 9Q 2 2Qii22 + 9Q11Q1111 - I8Q12Q1112 + 9Q11Q1122 - I8Q12Q1222 - 9Q22Q1111 
— 3Qinii2 + 6Q111222 + 9Q122222 — 9Q22Q1112 + 9Q11Q1222 — 9Q22Q1222 — I8Q12Q1122 + 9Q11Q1112 — I8Q12Q2222 
— Q111111 + 6Q111122 — 9Q112222 + 9Q 2 2Qn22 + 3Qn(5iiii — 6Q12Q1112 — 9Q11Q1122 + I8Q12Q1222 — 3Q22Q1111 



— 3Q111112 
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+ 3Q22Q1222 + 6Q12Q2222 
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